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Abstract 

We increase the efficiency of a recently proposed time integration scheme for time 
dependent quantum transport by using the Lanczos method for time evolution. We 
illustrate our modified scheme in terms of a simple one dimensional model. Our 
results show that the Lanczos-adapted scheme gives a large increase in numerical 
efficiency, and is an advantageous route for numerical time integration in ab initio 
treatment of open boundary quantum transport phenomena. 
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1 Introduction 

In many physical phenomena, practical limitations hinder a complete knowl- 
edge of all the degrees of freedom involved. Nanoscience has adopted such 
apparent shortcoming as its central paradigm, by exploiting the notion of a 
small system coupled to a macroscopic environment. A case in point is repre- 
sented by quantum transport phenomena, where two (or more) macroscopic 
leads are connected to a small central device (quantum constriction). 
Theoretical approaches to quantum transport can be broadly grouped in two 
categories, those based on a steady state formulation and those using a time 
dependent framework. Another discerning criterion can be the type of method 
used. In this case one can primarily distinguish among ab initio or model 
Hamiltonian methods. Finally, one can also consider a distinction based on 
the mathematical technique used: nonequilibrium-propagator, linear-response, 
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wavefunction-scattering, etc. Here we consider the time dependent quantum 
transport (TDQT) approach, which permits to follow the system during its 
time evolution after a bias has been applied. In this way, steady-state, tran- 
sients and a.c. currents can all be considered on equal footing and, in the pres- 
ence of dissipation, history dependence (memory effects) are also accounted 
for. An early formulation along these lines was introduced almost three decades 
ago [1]. 

For a quantitative description of TDQT, as for example required to obtain 
a theoretical figure of merit of the transient response of a real device, a de- 
scription at the ab initio level is certainly required. For this, one can resort 
to Time Dependent Density Functional Theory (TDDFT) [2,3]. In TDDFT, 
the TDQT problem is rigorously mapped onto a fictitious independent par- 
ticle problem. A formulation of TDQT within TDDFT has been introduced 
recently [4] . The practical applicability of the method has also been shown [5] , 
and the formulation has been extended to include classical nuclear degrees of 
freedom [6]. 

The purpose of this short communication is to show how the Lanczos algorithm 
for time evolution [7] can be applied to the case of open geometries as those 
encountered in time dependent quantum transport. This is done introducing a 
modification to the approach given in Ref. [5]. After a quick presentation of the 
Lanczos algorithm, we will review the method in Ref. [5]. Then we present our 
Lanczos-adapted method, and show comparative results for a model system, 
followed by some conclusive remarks. 



2 The Lanczos method 

We briefly summarize the Lanczos method, as given in [7]. A useful compar- 
ative study between the Lanczos method and other integration schemes can 
be found in [8]. Consider a system described by a TD Hamiltonian H(t). If, 
for example, we use the mid-point approximation for the time propagator and 
wish to evolve the system in the time interval (t + A, t), we obtain 

|<W = e-^ t+A ^ A \^ t ) (1) 

where |$ t ) is the (known) initial wavefunction. Consider a finite Lanczos se- 
quence {|Vfc)}, obtained by starting acting on the 'seed' = \V ). Using 
as a truncated basis, we get 

|^+a) « E l^> (Vk\e- lHLt \V ), (2) 

fc=0 

where Hz, is the tridiagonal representation for H(t + A/2) in such a basis. 
Inserting a complete set of eigenstates for the truncated space, H L \\) = e\\X), 
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where I^+a) is finally expressed in the basis of the original many body Hamil- 
tonian. The method requires a partial orthogonalization on the fly of the 
Lanczos basis in order to preserve accuracy along the trajectory. For a simple 
estimate of the truncation error in Eq.(2), see the discussion in [7]. 



3 Boundary Conditions in Time Dependent Quantum Transport 



An effective and viable strategy to TDQT is to consider large but finite sys- 
tems. Via an initial charge imbalance, a quasi-steady state current can be 
established, as clearly shown either in presence of electron-nuclear interac- 
tions [9] or when only electrons are considered [10]. A different approach, the 
one we consider here, is based on an open boundary formulation of the prob- 
lem [1,4], with a central region connected to two semi-infinite leads [5]. This 
approach has also been used in a mixed quantum-classical scheme to deal 
with electron -phonon systems in quantum transport, where the phonons are 
treated as classical fields (Ehrenfest Dynamics, ED) [6]. It is for this latter ap- 
proach, which has recently received some attention in the literature, that we 
present a Lanczos-adapted numerical scheme. We will consider for simplicity 
the purely electronic case: classical nuclear degrees of freedom can be added in 
a straightforward manner. Finally, a time-dependent embedding scheme has 
been considered very recently also in [11]. 



4 Time Evolution for Quantum Transport 



We provide here a brief presentation of the open boundary algorithm of Ref . [5] . 
No attempt of completeness is made and we refer to the original paper for a 
detailed derivation. In the following, we present the main formulas in the 
case of a strictly ID system (i.e. the leads have no translational invariance 
in the transverse direction), to provide the background needed to introduce 
our Lanczos-adapted scheme. The Hamiltonian we consider is H tot (t) = H e ; + 
W(t), where W(t) is the external perturbation. In a TDDFT approach, the 
initial, ground state is a single Slater determinant \^ g ). It is useful to divide 
the (ID) space into three regions. With s a the site label, we have the region 
L (corresponding to the left lead, with s < — (M + 2) ), the central region C ( 
with \s\ < M + 1 , i.e the device region contains 2M + 3 sites), and the region 
R (corresponding to right lead, with s > (M + 2)). The general structure 
of any bound, extended or resonant one particle eigenstate tj) in the Slater 
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determinant \^f g ) can be written as 



l +e -ikis + i_ e ikis s < _ M _ 2 

ip(s) \s\ < M+ 1 , 

R + e ikrS + R_e~ lkrS s > M + 2 



(4) 



To describe quantum transport, one needs to evolve in time the ground state 
configuration \^f g ), i.e. each one of the single particle eigenstates if) above. 
Introducing the projection operators Pl,c.r (for example, Pl = J2s&l I s ) 
we can write (/3 = L, C, R), for the generic single particle state, 



/3 



1^) 



(5) 



In the same way, we can project the Hamiltonian in the different regions 



H = ^H^, 

PP' 



H/3/3' = P«HPfl'. 



(6) 



Separating the contribution from the leads in W, the set of one-particle equa- 
tions becomes 



i~\iP(t)) = [n(t)+W leads (t)] |V(t)>, 



(7) 



with H(t) = H e i + W cc -(t), where H el is the electron one particle Hamilto- 
nian and Wc^(t) is the external potential projected in the central region C. 
Assuming metallic electrodes [5], 
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8 s , s >W R (t) s>M + 2 
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In the numerical time propagation, the time is discretized: t m = 2m5 , where 
5 is the timestep, m is an integer, and the explicit prefactor 2 is introduced for 
convenience in the formulas. In [5], the one-particle eigenstates are propagated 
from t m to t m+ i using a generalised Crank-Nicholson scheme. For the time 
evolution of each one of the one-particle states in \^ g ), one gets [5] 
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where |^ m ) = \%fj{t m )) and 
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H< m > = H cl + i [W cc (t m+1 ) + W cc (t m )} (10) 
wL m i = \[W leads (t m+1 ) + W leads (t m )}. (11) 



4-1 Propagation of One-Particle Eigenstates 



Using Eqs.(5,6), and after some algebra, the closed equation for the time- 
evolution in the central region is 



l4" +1) > = \ c - Zlt) \^ ] ) -m E % (l7?>> + ic£°>) , (12) 



where 



win) _ (13 ) 
Q "1 + ^' (13) 
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nf^IM 2 . ( 14 ) 
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and 



H$=H<&-i* £ H c 1 K aC = lI&-i6 £ B?. (15) 

The B^ ** matrices have only one non-zero element, 

B (o)] s ^ fe (o)|VM-i^-M- a « = L^ 

[ ^s,Af+l^s,M+l a = R 

with = ~ 1+v/ 2 1 <52 4 ' 5 ~^~ an< ^ ^ ^he hopping parameter in the leads. The ex- 
pression for the source state \^') and the memory state |Ci™' ) ) are [5]: 



\i { : ) )=g { : ) - — l —^K) (is) 
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where \u a ) is a unit vector such that 



(s\u a ) = { . (19) 

8 s ,M+i a = R 



The scalar quantities Zf™ and a = L, R are given by 



^ + ^^'^O 0«l^' +1) > + («a|^>) , (20) 



2* £S ^ 

v y (1 — 2i<5 cos (z a )) 

+ (a + e iz ^ M+ V + a„e-^ M+1 A x idT ^ ~ — — ■ + 

j=o (1 — 2i<5cos (<2 a )) n J v 

and Zq, = k[ for a = L while z a = fc r for a = R. For n > 2, the quantities fr*-™* 1 
in the Eqs. (20,21) are obtained by recursion: 



= _ 5 2 -^ C221 

n-l (h(i) + _L ftO'- 2 )^ ft(«-2-i) 
-£ 2 V A ^ 

^ 1 + 25 2 6(°) 

and 6( n< °) = 0,feW = }; 2 gg feW and &(°> the same as in Eq.(16). 



5 Lanczos-adapted algorithm 



The basic idea behind the algorithm illustrated in the previous Section is to 
discretize the time axis via the Crank- Nicholson algorithm before perform- 
ing the partitioning in L, C, R regions [5]. One could devise doing the same 
for the Lanczos algorithm; however, noncommuting parts of the Hamiltonian 
would appear in the exponent this time, rendering formal manipulations more 
involved. Here, we consider a simple shortcut that, while improving the nu- 
merical efficiency of the algorithm in [5], has the same degree of accuracy( i.e. 
it is second order in 5) but avoids working with the Lanczos scheme before 
the partitioning. Looking at Eq.(12), we notice that the explicit action of 
occurs in two specific terms: 
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1*)= " ( 23 ) 
IX2> = -^o l«a) (24) 

where is the contribution to IV'c ) f rom the central region, and \x2) 
enters the expressions for the source and memory states. For since 5 — > 0, 
one can write, up to order two in S 



IXi>= ^"l°g ) l^ ) >«e-^|^ ) >. (25) 



For the case of IX2), we define the following quantities: 



A ± = (26) 



which permit to rewrite [xz) as 



u a ) + 0(5 3 ) (27) 

If necessary, one can go to higher orders, by imposing that (1 + <5x) _1 = 
A + J2k e akSx and finding the coefficients A, {ctk} by comparison of the two 
expressions order by order in S (in general, the {0^} will be complex). We 
note that the same Lanczos sequence of basis vectors is required for both 
exponentials in Eq.(27). 

All terms which appear in the propagation scheme [5] and that involve H^j , 
have been re-expressed in terms of exponentials, so that Lanczos propagation 
can be used; finally, since H^j is complex, Eq.(15), it is convenient to split 
the exponentials; for small 5, 
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e -w=5 , « e -*'E.^ e -^ e -' a E.^ (28) 

e -*A±H$ ^ e -|A ± E Q BL 0) e -iA ±H W e -f A ± E a BL°» 

(29) 

For the ID case, the advantage is immediate: the B^ * 1 in Eq.(16) have only 
one non- vanishing entry and the outer exponentials in Eq. (28,29) reduce to 
scalars (here, we do not address the 3D case; however, we expect that the 
splitting will still provide a simplification). 
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Fig. 1. Results for the current for different numbers of sites in the central region, 
for a bias Ul = 0.5|V|. Color coding is specified in the top right inset panel. 



6 A Numerical Example 



The advantage of the modified scheme presented here is the possibility of ma- 
nipulating very efficiently exponentials (via the Lanczos scheme), thus being 
able to deal with larger scale problems in a faster way. We have performed 
some tests for a simple spinless model of QT, namely a 3D central region 
C connected to two semi-infinite ID metallic leads, in the half-filling regime. 
The leads are described by a nearest neighbour, tight binding Hamiltonian 
(nnTBH) with hopping parameter V = — 1. The central region, as shown in 
Fig.l, is made of a short central chain of five sites connected with two identi- 
cal clusters. Such clusters are composed by periodically repeated layers, each 
layer containing four atoms arranged in a square. For technical reasons, the 
rightmost (leftmost) of the left (right) lead is also included in C. The single 
particle Hamiltonian in C is also a nnTBH, where Vc = —0.3. The number 
of Nc sites in the central region is Nq = + 7, and we vary it by chang- 
ing Nl. If we increase Nl, we can think of our system as a five-site chain 
connected to finite, but progressively longer 3D leads (the latter are in turn 
connected to the ID, truly semi-infinite leads). In our QT simulations, there 
will be a transient, but increasingly longer, time interval before the truly ID 
nature of the electron reservoirs will manifest. We have analysed the current 
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Table 1 

Comparison between different schemes of numerical integration for the algorithm of 
Ref. [5] for central regions of different size. Execution times are in arbitrary units. 

N L 5 10 15 20 30 

Inversion 1.00 1.51 2.14 3.30 14.94 
Lanczos 1.55 2.02 2.50 2.99 4.02 

j = 2V J2° k cc Im[ip k (R)iljl(R + 1)] at the central site R = (the grey-shaded 
circle in Fig.l) as a function of Nl- For any fixed time i, on increasing Nl the 
current j converges to a specific value; deviations from the converged value 
occur at longer times for greater values of Nl, because the ID nature of the 
real reservoirs enters at later stages for longer clusters. To assess the efficiency 
of our modified scheme, we have calculated the currents of Fig.l for different 
Nl in two ways (which differ on how Eqs.(23, 24) are computed). Namely, we 
i) used standard LAPACK routines to compute the inverse of the operator 
lc + idHtff , and ii ) used the Lanczos-adapted scheme introduced here. We 
note that, in analogy to [12], another way to manipulate Eq.(23) is to iii) solve 
a linear system, after recasting Eq.(23) as 



(lo + i*Hff|xi> = (lc-im^)\^). (30) 

Such linear system is to be solved for each single particle state, and this is 
expected to become computationally unfavourable (unless the Hamiltonian 
has a special structure such as band-diagonal, sparse, etc.) when the number 
of single particle states in the Slater determinant and/or the size of central 
region become large. On the other hand, the operator in Eq.(23) is state- 
independent, and the inversion can be performed before entering the loop 
for the single particle states in the Slater determinant. Accordingly, we did 
not consider iii) in our numerical comparisons. In all calculations we used a 
timestep 5 = 0.01 1 ~ x , with Ai=5000 timesteps. For the short iterated Lanc- 
zos scheme, we used 6 iterations/timestep. Results for the execution times, 
as a function of Nl are shown in Table 1. We see that on increasing Nl, the 
Lanczos adapted scheme becomes significantly more efficient than i). We ex- 
pect this to be a general trend: for genuine 3D systems/leads, the advantage of 
a Lanczos-adapted time evolution should then become even more significant. 
At the same time, the actual figures of relative numerical efficiency between i) 
and ii) in Table 1 should be considered only as indicative, since we have not 
performed a careful optimization of the Lanczos-adapted algorithm/code (an 
optimized code could further improve the numerical performance). 
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7 Conclusions 



In this short note, we have described a simple way to increase the numeri- 
cal efficiency of a recently proposed algorithm for time dependent quantum 
transport. We tested the efficiency of the proposed scheme in terms of a model 
system. While our modifications to the original algorithm are rather simple, 
we expect that the practical advantage of such modifications to be significant, 
since future time dependent ab initio calculations for quantum transport in 
realistic structures are expected to involve sizeable active regions, i.e. large 
configuration spaces and large scale calculations. This work was supported by 
EU 6th framework Network of Excellence NANOQUANTA (NMP4-CT-2004- 
500198). 
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